Efficient implementation of a van der Waals density functional: Application to 

double-wall carbon nanotubes 
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We present an efficient impiementation of the van der Waais density functional of Dion et al 
[Pfiys. Rev. Lett. 92, 246401 (2004)], wliidi expresses tfie nonlocal correlation energy as a double 
spacial integral. We factorize the integration kernel and use fast Fourier transforms to evaluate 
the selfconsistent potential, total energy, and atomic forces, in log A'' operations. The resulting 
overhead in total computational cost, over semilocal functionals, is very moderate for medium and 
large systems. We apply the method to calculate the binding energies and the barriers for relative 
translation and rotation in double-wall carbon nanotubes. 

PACS numbers: 31.15.eg, 71.15.-m, 61.46.Fg 



C/2 



X3 

o 



> 

O 

OO 

o 



X 



Density functional theory (DFT) has become the 
method of choice for first-principles simulations of static 
and dynamical properties of complex materials with 
strong ionic, covalent, and metallic interactions. How- 
ever, weak van der Waals (vdW) interactions are also 
essential for many systems and processes, like molecular 
solids and liquids, surface adsorption, and biological reac- 
tions [l| . Local or semilocal density functionals obviously 
cannot describe asymptotically the nonlocal dispersion 
correlations. At binding distances, they have been fre- 
quently found to give reasonable results 0, Q but their 
ability to do so is generally very sensitive to the spe- 
cific functional used and its parametrization details, what 
makes the "a6 initid^ character of this approach rather 
questionable. Thus, the simulation of vdW systems has 
typically relied on atom-atom potentials with the conven- 
tional [J| asymptotic behavior and with parameters 
fitted to empirical data or to accurate quantum chem- 
istry calculations of simple molecules. Such potentials 
are also added as plug-ins to ab initio semilocal density 
functionals [1, @|. Another approach includes vdW in- 
teractions through effective atom-electron pseudopoten- 
tials 0]. However, the accuracy and reliability of such 
approaches is limited because vdW energies arise from 
electron-electron correlations that depend not only on the 
atomic species but also on their chemical environment. 
More ah initio wavefunction-dependent approaches are 
more reliable but also much more expensive Q . 

Thus, a key development has been the proposal by 
Dion et Q of a universal nonlocal energy functional 
of the electron density n(r) with the form 



E^M^)] = E^^^[n{v)] + E^^^Hv)] + E:^[n{r)] (1) 

where the exchange energy E^'^^ is described through 
the semilocal generalized gradient approximation 
(GGA) and the correlation energy has a local part 
E^^^, described in the local density approximation 



(LDA), and a nonlocal (nl) part E" given by 



E. 



nl 



Hr)] = ^ 



(fri d?V2 n{vi) n{r2) 0(51,92,^12) 

(2) 

where ri2 = |ri — r2|, and gi, (72 are the values of a uni- 
versal function qo{n{r), |Vn(r)|), evaluated at ri and r-z- 
The kernel has also a precise and universal form that 
in fact depends only on two variables di ~ qiTi2 and 
d2 = 92^12, but it can obviously be written also as a 
function of Qi, 52, and ri2, what we will find convenient. 
The shape of obeys that: i) i?"' is strictly zero for 
any system with constant density; and ii) the interac- 
tion between any two molecules has the correct de- 
pendence for large separations r. Using a direct eval- 
uation of Eq.(l2]), this vdW functional has been applied 
successfully to a variety of systems, including interac- 
tions between pairs of atoms and molecules, molecules 
adsorbed on surfaces, niolecular solids, and biological sys- 
tems. i,[ii|,|Tl,[ii,|il 

If qi and q2 in Eq. ([2|) were fixed values, independent 
of ri and r2, i?"' would be a simple convolution, like 
the Coulomb energy, that could be evaluated by Fourier 
methods. Therefore, our key step for an efficient imple- 
mentation is to expand the kernel (j) as 

0(91, 92, '"12) ~ ^(^(9a,9/3,fl2) Pa{qi) Pp{q2) (3) 

where qa are fixed values, chosen to ensure a good in- 
terpolation of function cf). In order to illustrate how the 
factorization ([3]) can be performed in a systematic way, 
we consider first the interpolation of a function f{x) us- 
ing a linear scheme, like those of Lagrange, Fourier, or 
splines: 



f{x) ^^fa Pa{x) 



(4) 



where fa = f{xa) and Pa{x) is the function resulting 
from the interpolation of the particular values fp = 



2 



5ai3- In Lagrange interpolation, it is a polynomial of 
given order. In Fourier interpolation it has the form 
sin(7r (x— ) / Ax) / (7r(x— Xq, ) / Ax) . We use cubic splines, 
in which Pa{x) is a succession of cubic polynomials in ev- 
ery interval [x/3,a;^+i], matching in value and first two 
derivatives at every point xp. Notice that Pa{x) depends 
on the interpolation scheme and on the (fixed) points Xa , 
but not on the interpolated function. In two-dimensional 
interpolation, one typically interpolates first in one vari- 
able and then in the other: 



f{x,y) ~ ^f{x,yp) pp{y) 



- ^y^fi^a^yp) Pc.{x)\pi3{y) (5) 

what shows that such an interpolation leads automati- 
cally to an expansion in terms of factored functions of x 
and y. Thus, Eq. ([3]) is just the interpolation of a three- 
dimensional function in its first two variables. In this 
latter case, however, the interpolation points qa must be 
appropriate for every value of the third variable ri2. 

The fact that ri2 acts as a scaling factor (i. e. increas- 
ing ri2 merely "contracts" (/) as a function of qi and (72, 
without changing its shape) suggests a logarithmic mesh 
of points qa, in which [qa+i ~ la) = ^{(la — qa-i), with 
A > 1. Such a logarithmic mesh is also suggested by the 
shape of 4>{di,d2) shown in Fig. 1 of ref. [91]. We have 
found that Na ~ 20 interpolation points qa are sufficient 
for an accurate description of </) up to a cutoff qc at which 
we artificially "saturate" the original function (70(^1 l^"-l) 
by redefining 

C*(n,|Vn|) = /i[qo(n,|Vn|),g,] (6) 

where h{x, Xc) is a smooth function such that h{x, Xc) — 
X for X < Xc and h(x, Xc) Xc for x 00: 



h{x, Xc 



Xc 



1 — cxp 



jx/XcY' 

■m—1 



m 



(7) 



with TUc ~ 12 and qc ^ 5 a.u. Higher qo values arc ob- 
tained only for very large n(r) (i.e. close to the nucleus, 
where i?"' is negligible compared to other terms in Exc), 
and for large |Vri|/n (in the electron density tails, where 
i?"' is negligible because of the factor n{r) in the inte- 
grand of Eq. ([2])). In what follows, we will omit, but 
assume, supcrindcx "sat" in qa{n, |Vri|). 

A minor but significant difficulty is that (j){di,d2) has 
a logarithmic divergence when di,d2 — > 0, what prevents 
its straightforward interpolation. Therefore, we interpo- 
late and use instead a modified "soft" form 



4>s (^1,^2) 



9o + <P2C 
4>{di,d2) 



(pid'^ if d < ds 
otherwise. 



(8) 



where d = y/d^+d!^. 0o a-nd ds are fixed parameters, and 
(/)2, 04 are adjusted so that (jysidi, ^2) and 0(di, ^2) match 



in value and slope at d = ds (for given ^2/^1). This 
modification leads to a change in i?"' , which is corrected 
using a local density approximation: 



ae:: 



(9) 



where 



AeL''(r) = ^ / A7Tr'^dr'[cl,{q,q,r')-Mq,q,r')] 



n{r) 
'2^ 



And^dd [(f>{d, d) - (j>s{d, d)] (10) 



with q = qo{n{r), Vn(r)). The evaluation of AE"^ and its 
derivatives is performed, like that of the semilocal terms 
in Eq. ([1]), as in ref. jT^. In what follows, we will assume, 
but omit for simplicity, the subindex s in 0^. 
Substitution of jS]) into ([2]) leads to 



^ [ I ^'^1 ^^'^2 0a(ri) ep{r2) Mri2) 

= IT./"^'^ ^/3(k) q^Mk) (11) 



a/3 



where 0a(r) = n{r)pa{qo{n{r),Vn{r))) and 0a(k) is its 
Fourier transform. Equally, (papik) is the Fourier trans- 
form of (jiapif) = 4'{qa,q0,T'). It can be calculated in 
spherical coordinates, and stored in a radial mesh of 
points k for convenient interpolation. Thus, the heavier 
part of the calculation is the fast Fourier transforms of 
the Na functions Oa (r) , which still have a very moderate 
cost in a typical density functional calculation. 

The evaluation of atomic forces requires the use of the 
Hellman-Feynman theorem, which holds only if the full 
energy functional is minimized selfconsistcntly. In turn, 
this requires the nonlocal part of the correlation poten- 
tial, i. e. the functional derivative of Eq. ([2]) [l6|. To 
handle the gradient dependence in go("-5 Vn) we use the 
same technique as in ref. [l5| : approximating the spatial 
integrals by sums in a uniform grid of points, and the gra- 
dients by finite differences in the same grid. This makes 
i?"' an ordinary function of the densities Ui at fixed grid 
points Ti , allowing to perform conventional partial deriva- 
tives, rather than functional derivatives. Besides its con- 
ceptual simplicity, this method ensures a perfect consis- 
tency between the calculated potential and the energy: 

E^' = iAr!^ ^ Cl>ap{n,) (12) 

ap ij 

where AD, is the volume per grid point and 9ai = 
niPaiqaini^Vrii)). Notice that (j^apifij) does not depend 
on rii, since the values qa are fixed. A straightforward 



3 



derivation then gives 
1 dEV-^ 



0.1 r 



Ail dui 



E 



dn. 



(13) 



where dS/njjdni are fixed coefBcicnts (determined by the 
finite difference formuia used for VtIj) that depend only 
on Vij and that are nonzero only for small r.y . Also, 

Uai = Aft^^^^Opj (paliirij) (14) 
/3 3 

is a convolution that can be obtained using fast Fourier 
transforms since (apart from tt and volume factors) 



d^V2 0l3{y^2) (f>a0iri2) 



d-'k e"'--^ e^ik) 0„^(fc). 



(15) 

Thus, a selfconsistency step requires direct trans- 
forms to find 0a{k) and inverse transforms to ob- 
tain Ua{r). The calculation of the atomic forces does not 
require any additional effort, since the nonlocal contribu- 
tion w"' is simply added to the semilocal terms [isj in vf^ 
and to the rest of the effective potential. Notice that the 
implementation is independent of the basis set, accept- 
ing n{ri) in a uniform real space grid and returning 
Exc and Vxc{j^i) in the same grid. It has been checked 
that it reproduces accurately the results obtained by di- 
rect evaluation of Eq. ([2]) (and, eventually, its functional 
derivative [31) for a variety of systems (ol. [l^. 

We have applied the above method to study the inter- 
action between the concentric layers of double-wall car- 
bon nanotubes (DWNT). Such interactions are crucial for 
different nanodevices proposed recently 17 , T3| and they 
have been studied with semiempirical potentials [Toj and 
with a local DFT functional [20, Hll, [2§| . We have used 
the SIESTA code [H, Q with an optimized ^ triple- 
(■-l-polaratization basis set of pseudoatomic orbitals, cor- 
recting for basis set superposition errors (BSSE). The 
integration grids in real and reciprocal space had cutoffs 
of 300 Ry and 20 A, respectively. The cutoff parameter 
for fc-point sampling [2^ was 20 A, ensuring at least 34, 
20 and 14 fc-points for the armchair, zigzag and chiral 
DWNTs studied, respectively. The atomic forces were 
relaxed to less than 20 meV/A. 

Figure [1] shows the calculated interaction energy be- 
tween two rigid SWNTs, relaxed independently, as a 
function of their interwall separation (difference of radii) . 
It also shows the DWNT formation energies, defined as 
the difference between the total energy of the relaxed 
DWNT and that of the two SWNTs. The calculated 
tubes (m,m)@(n,n) (armchair), (m,0)@(n,0) (zigzag) and 
(8,2)@(16,4) (chiral) were chosen for their conmensura- 
bility in the longitudinal direction, as well as for com- 
parison with prior calculations. It can be seen that the 
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FIG. 1; Interaction and formation energies between differ- 
ent double wall carbon nanotubes, as a function of inter- 
wall separation, using LDA [23] (squares), GGA [2^ (circles), 
and van der Waals (triangles) functionals. Interaction 
energies (empty symbols) are between two individually re- 
laxed rigid tubes. Formation energies (filled symbols) include 
also the geometry relaxation induced by the interaction, that 
modifies their interwall separation. The tube geometries are 
(5,5)@(n,n) (0,0, A)), (m,m)@(n,n) m>5 (v), (m,0)@(n,0) 
(>), and (8,2)(Q)(16,4) (<). For comparison, we also show the 
interaction energies for two fiat graphene layers (lines). All 
energies are divided by the total number of atoms in both 
tubes. 



interaction energy depends neglegibly on chirality and 
curvature, being very well represented by the interaction 
between two flat graphene layers. On the other hand, 
the formation energy, that includes the relaxation of the 
radii induced by the interaction, shows a steeper repul- 
sion than between flat graphene layers. In agreement 
with previous results for graphene and graphite, we find 
that the LDA works reasonably well. For sufficiently long 
tubes, in which the border effects can be neglected, the 
calculated vdW interaction energy gives a telescopic con- 
traction force [3| F = 0.91N/m x d, where d is the mean 
of the inner and outer tube diameters. 

Next, the two concentric tubes of the DWNT were 
moved rigidly, relative to each other, in order to con- 
struct rotation-translation energy maps. To generate 
these maps, we first project the inner tube coordinates 
onto the outer tube surface, i. e. we multiply its x and y 
coordinates (the tube axis being z) by the ratio Rout /Rm 
between the two radii. We then unroll the coordinates 
of both tubes onto a fiat surface, repeating them peri- 
odically also in the x axis. This gives two flat periodic 
lattices (conmensurate in the cases considerd) with re- 
ciprocal unit cell vectors and hi,i = 1,2. The energy 
maps can then be represented, as a function of the po- 
sition X on this surface, relative to the minimum, by an 
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in generating basis sets and testing our implementation. 
This work has been founded by grant FIS2006-12117 from 
the Spanish Ministery of Science. 



TABLE I: Periodicities (Aii = 27r/Gi, in A) and energy 
barriers Ui (in meV per outer tube atom) for translation (i = 
z) and rotation (i — <p) of the outer tube, relative to the 
inner tube, in double wall carbon nanotubes. Ax^ lengths 
are along the outer tube circunference. For the (8,2)@(16,4) 
tube we found that all the LDA and vdW barriers are smaller 
than out computational accuracy of ~ 0.01 meV/atom. 



expansion of the form 

C/(x) = f/o - i ^ C/g cos(G • x) (16) 

where G are the superlattice wavevectors, common to 
the reciprocal lattices a and b, and Uc, are the barrier 
heights for motion along G. We have found that limit- 
ing this expansion to the first two wavevector stars, ±Gi 
and ±G2 (which, in the cases studied, are parallel and 
orthogonal to the axial direction), gives a good approx- 
imation to the cases studied, with the parameters given 
in Table U 

Overall, the relative values of these barrier heights are 
in qualitative agreement with previous calculations, i. 
e. larger barrier distances lead to larger barrier heights. 
Quantitatively, however, those calculations vary by an or- 
der of magnitude depending on the models used [3, [2^ . 
Our calculated LDA barriers are similar to those of 
refs. [1^, [m, [111. 



The small discrepancies with ref. [22| 
may be due to the different basis sets and to our finer 
fc-point sampling. Again, we find that the LDA does 
a rather good job for these systems, compared to the 
more complex vdW functional. Nevertheless, we find 
that LDA systematically underestimates the interaction 
energies and that it overestimates the barrier heights, 
relatively to the vdW results. 

In conclusion, we have described an efficient algorithm 
to include van de Waals interactions through the selfcon- 
sistcnt treatment of a nonlocal ab initio functional pro- 
posed recently. Q Typical overheads in the total compu- 
tation time, using the SIESTA code, over that required 
by LDA or GGA, are a factor ~ 5 in a two-atom system 
and a ~ 10% increase for ~ 150 atoms. Using this imple- 
mentation, we have calculated the interaction energies, 
as well as the barriers for relative displacement, between 
concentric tubes in several armchair, zigzag, and chiral 
DWNTs. 

We would like to thank E. Anglada, E. Artacho, and 
D. C. Langreth for many discussions and for their help 
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